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ABSTRACT 

Large  Eddy  Simulation  (LES)  is  developed  in  a  general  purpose  Finite  Element  code.  It  is  first  shown  that  the  standard 
numerical  schemes  developed  for  the  Reynolds  Averaged  Navier  Stokes  Equations  (RANSE)  are  not  suitable  for  LES. 
Validation  of  the  final  centred  and  collocated  scheme  is  carried  out  for  grid  turbulence  and  channel  flow.  The  same  code 
is  then  applied  to  a  complex  flow  that  has  challenged  RANSE  modelling  for  years:  a  crossflow  in  a  tube  bundle.  In 
opposition  to  such  models,  the  LES  yields  results  in  good  agreement  with  the  experimental  data  of  Simonin  et  al., 
including  for  the  Reynolds  stresses. 


1.  INTRODUCTION 

Large  Eddy  Simulation  (LES)  has  long  been  restricted  to 
academic  flows  (homogeneous  or  channel  flows),  while 
models  for  Reynolds  Averaged  Navier  Stokes  Equations 
(RANSE)  are  commonly  applied  to  three-dimensional 
flows  in  complex  geometries.  Some  fairly  complex 
domains  are  simulated  by  LES  :  flow  around  a  cube  or  a 
square  cylinder  (1).  However  LES  simulations  require 
high  mesh  refinements  near  the  body  and  induces  in 
structured  meshes  a  very  high  number  of  mesh  points. 
The  advantage  of  unstructured  grids  is  to  allow  local 
mesh  refinement  (2).  But  on  the  other  hand,  the  highly 
accurate  schemes  used  for  LES  (Pseudo-spectral  or  Fade 
schemes)  are  applicable  only  to  cartesian  grids.  Only 
two  applications  of  LES  code  on  unstructured  grids  for 
incompressible  fluid  developments  are  known  (2)  &  (3). 
We  present  here  the  development  of  a  LES  code  for 
incompressible  fluid  through  a  Finite  Element  Method 
(4).  This  LES  version  is  based  on  the  N3S  code 
developed  by  Chabard  et  aL{5). 

The  organisation  of  this  paper  is  as  follows.  The 
numerical  issues  which  were  explored  are  described  in 
Sec.  2.  The  effect  of  numerical  method  on  LES  of 
channel  flow  is  analysed  in  Sec.  3.  The  final  numerical 
tool  is  applied  to  tube  bundle  flow  in  Sec.  4. 
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The  filtered  Navier-Stokes  equations: 
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The  Smagorinsky  model  is  a  simple  mixing  length 
model: 
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with  Vf  the  turbulent  viscosity 

and  A  the  characteristic  length  scale,  A  =  2/2 , 
with  h  the  mesh  spacing  for  a  regular  mesh 
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C=Cf ,  is  the  Smagorinsky  constant 


The  subgrid-modelling  is  an  important  problem  for  LES, 
but  we  focus  here  on  the  numerical  requirements  of  LES 
on  unstructured  grid. 

2.1.  Time  advancement 


2.  NUMERICAL  ISSUES 

In  opposition  to  the  Reynolds  Averaged  Navier-Stokes 
Equations  (RANSE),  with  LES  one  explicitly  computes 
the  large  scales  of  motion  and  introduces  a  model  for 
only  the  small  motions.  LES  is  based  upon  the 
application  of  a  spatial  filtering  operation  (denoted  by  a 
bar)  to  the  three-dimensional  unsteady  Navier-Stokes 
equations.  We  assume,  and  subsequently  check  from 
energy  spectra,  that  the  filtering  is  induced  by  the 
numerical  scheme  and  the  computational  grid. 

Let  u  the  (filtered)  velocity  of  large  scale  motion.  The 
subgrid-scale  stresses  (SGS)  result  from  the 

influence  of  unresolved  scales  on  large  scales,  and  are 
defined  by  Ty  =«,  Uj  -  Uj 


Two  different  methods  have  been  tested:  the  method  of 
characteristics,  and  the  centred  Adams-Bashford,  Crank- 
Nicholson  scheme  (AB-CN). 

The  characteristics  scheme  deals  only  with  the  advection 
step^and  is  used  within  a  fractional  time-step  method. 
For  pure  advection,  velocities  are  conserved  along  the 

dx 

characteristic  curve  defined  by:  — =m, 

dt 

=  v  —  rrw.  At  each  time  step  t",  this  curve  is 
dt  dt 

followed  backwards  from  each  node.  Once  the  location 
of  the  fluid  particle  at  t"  ‘  is  known,  the  velocity  is 
interpolated  by  quadratic  or  higher  order  interpolation 
functions.  Then,  the  Stokes  problem  will  be  solved  for 
this  intermediate  velocity  field. 


This  approach  is  very  attractive,  because  it  introduces  a 
natural  upwinding  and  is  unconditionally  stable:  it  is 
used  daily  for  k-£  calculations.  For  steady  state 
problems,  it  has  very  low  numerical  diffusion  since 
cubic  or  hermitian  interpolation  can  be  used  at  the  end 
of  the  characteristic  curve.  But,  for  a  highly  unsteady 
flow  as  in  LES,  and  since  the  discretised  characteristics 
method  mixes  interpolation  in  space  and  time,  the  first 
order  temporal  interpolation  results  in  unacceptable 
numerical  damping.  As  for  any  upwind  scheme,  an 
artificial  viscosity  can  be  exhibited.  Namely,  for  the 
characteristics  method,  in  ID,  for  a  regular  grid  and  u 
the  advection  velocity: 

9/  A,-/" 

dt  At  At 

By  a  Taylor  expansion: 

3/  A  ,  u 

If 

At  =-j  lie  Ax ,  c  is  the  Courant  number. 

The  scheme  is  first-order  accurate  in  time  and 
equivalently  in  space  if  the  Courant  number  is 
prescribed,  this  whatever  the  order  of  the  spatial 
interpolation  at  the  end  of  the  characteristic.  This  was 
identified  in  (6)  and  resolved  for  cartesian  grids  by 
introducing  a  moving  frame  of  reference  at  each  time 
step,  but  cannot  be  extended  to  unstructured  grids. 

The  very  classical  Adams-Bashforth  Crank-Nicholson 
(AB-CN)  scheme  was  then  chosen  to  evaluate  the 
benefits  of  a  centred  scheme.  It  is  translated  by  the 
variational  formulation  to: 

Jp^  + jVu.  Vv^Q  = 

a  n 

-jfVxdQ.  -d-Jfv-i- V, jVu”. Vv^Q 

o.  Q 

“tC*”  -i-boundary  terms 

with  U  the  velocity  field  after  advection-diffusion  step 
and  V  the  velocity  basis  function. 

The  influence  of  the  formulation  of  convective  term  has 
been  investigates  in  detail,  but  not  in  the  FEM  context. 
Here  we  chose  the  straightforward  advective  formu¬ 
lation  (6):  C”  =  f  III  —^ydO. 

h 

Linear  basis  function  are  used  for  the  variational 
formulation:  the  scheme  is  fourth-order  accurate  for  the 
convective  term  and  second-order  for  the  diffusive  term 
(8).  However,  to  avoid  the  extraneous  computational 
cost  implied  by  the  mass-matrix,  this  matrix  is 
condensed  on  the  diagonal.  This  provides  a  second- 
order  accurate  scheme.  The  linear  analysis  proves  that 
this  scheme  is  stable  for  a  Courant  number  less  than  0.2 
These  schemes  were  compared  on  isotropic  turbulence 
decay.  Comparison  of  the  energy  spectra  with 
experiment  do  not  tell  much  since  the  dynamic  model 
automatically  decreases  the  Smagorinsky  constant  if  a 


dissipative  scheme  is  used  (but  this  is  only  acceptable 
for  grid  turbulence).  The  spectral  transfer  function  is 
more  significant.  The  variation  of  the  spectrum  of 
turbulent  kinetic  energy  during  the  advection  step  was 
evaluated  on  fifty  samples.  The  averaged  variation  is 
drawn  on  Fig.  1  for  two  meshes:  a  coarse  one  (N=32 
points)  and  a  refined  grid  (N=64  points).  The  spectral 
analysis  shows  that,  for  homogeneous  turbulent  flow,  the 
advection  implies  a  transfer  of  energy  from  the  larger 
structures  (small  wave-number  k)  to  the  smaller  ones 
(large  k).  The  numerical  test  shows  that  the 
characteristic  scheme  is  dissipative  on  the  mean  and 
very  dependent  of  the  grid:  the  scheme  is  first-order 
accurate  in  space  because  of  the  highly  unsteady  nature 
of  the  smaller  scales.  On  the  contrary,  the  centred 
scheme  effectively  transfers  energy  from  the  small 
wave-numbers  to  the  large  ones  and  fairly  independent 

of  the  spatial  resolution. 
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Fig.  1  Spectral  transfer  function  T(k)  in  grid  turbulence 

We  stress  here  that  the  paradoxical  choices  of  upwind 
biased  schemes  for  most  RANSE  computations  and 
centred  schemes  for  LES  should  be  analysed  in  the 
spatial-temporal  frame,  and  is  due  to  the  highly  unsteady 
nature  of  the  smallest  resolved  scales  in  LES. 

2.2.  Spatial  discretization 

Spatial  discretization  relies  on  tetrahedral  elements, 
which  enable  a  large  flexibility  in  meshing.  We 
restricted  ourselves  to  linear  shape  functions.  We 
introduce  here  the  PI -PI  element  based  on  collocated 
velocity  and  pressure  nodes  (fig.2)  (hardly  ever  used  in 
the  F.  E.  literature,  except  to  point  out  it’s  defects).  The 
standard  Finite  Element  discretization  (used  in  N3S 
code)  is  the  Pl-isoP2  element  (fig.2)  based  on  two 
levels  of  discretization:  the  pressure  is  interpolated 
linearly  on  the  macro-element,  while  the  velocity  is 
interpolated  linearly  on  each  sub-element. 

The  Pl-isoP2  element  verifies  the  well-known  "Tnf-sup" 
condition,  and  ensures  the  existence  and  uniqueness  of 
the  solution  of  the  discretized  Stokes  equation  (8).  On 
the  contrary,  the  PI -PI  element  does  not  verify  this 
condition  and  provides  spurious  pressure  modes. 


The  modes  are  damped  out  by  an  Arakawa  correction. 
These  modes  appear  mainly  when  the  grid  is  generated 
by  splitting  an  initially  cartesian  grid,  but  the  correction 
is  hardly  necessary  when  the  grid  is  fully  unstructured. 


Pl-isoP2  element  Pl-Pl  element 

velocity'  filter:  2h  v.  p.  filter:  2h 

pressure  filter:  4h 


Figure  2  Choice  of  Elements  (in  2D) 


Since  the  filter  is  induced  implicitly  by  the  grid  and 
numerical  scheme,  it  clearly  appears  that  the  filter-width 
is  twice  larger  for  pressure  than  for  velocity  in  the  case 
of  Pl-isoP2  element  (see  Fig.  2).  The  interactions 
between  resolved  and  subgrid-scales  are  modified  and  it 
may  be  necessary  to  develop  new  subgrid-models.  On 
the  contrary,  the  filter-width  are  identical  for  velocity 
and  pressure  in  the  case  of  Pl-Pl  element. 

The  treatment  of  incompressibility  is  essential  as  well. 
Brezzi  et  cz/.(9)  defined  a  constraint  ratio:  the  ratio 
between  pressure  degrees  of  freedom  and  velocity 
degrees  of  freedom.  In  3D,  this  ratio  is  equal  to  1/3  for 
Pl-Pl  element  and  to  1/4  for  Pl-isoP2  element.  Hence, 
the  Pl-isoP2  element  is  weakly  compressible.  Pelletier 
et  underlines  one  consequence  of  this  behaviour: 


llu-v,  <C 


"V,€V. 


^h^Qh 


with  u  and  p  the  solutions  of  the  continuous  Stokes 
problem,  and  p^  the  solutions  of  the  discrete 
problem.  The  precision  on  velocity  depends  both  on  the 
approximation  of  velocity  and  pressure. 

A  laminar  test-case  highlights  the  effect  of  the  low 
resolution  of  pressure.  In  the  ’Taylor-Green"  vortices 
flow,  the  shape  of  the  cellular  vortices  remain 
unchanged.  The  analytical  solution  is  known: 

U  --s\n{kx)  cos{ky)  e~^^^ 

<  V=  cos(/cjc)sin(/:y)^”^^^ 

P  =  j[cos(2/:jc)  +  cos(2^y)] 
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Fig.  3  Under-resolved  Taylor  vortices.  Blue  (light  grey): 
velocities.  Red  (dark  grey):  pressure  gradients. 


The  test  is  conducted  with  a  coarse  resolution:  8  points 
per  wave-scale  for  velocity.  After  50  time  steps,  the  L" 
error  on  the  velocity  component  u  is  0.08  for  Pl-Pl  and 
0.19  for  Pl-isoP2.  Figure  3  shows  that  the  angles  and 
vector  amplitudes  of  velocity  and  pressure  gradient  are 
well  captured  by  Pl-Pl  but  not  by  Pl-isoP2.  This  is 
very  important  because  the  energy  exchanges  between 
velocity  components  are  driven  by  velocity  pressure 
gradient  product,  as  in  the  Reynolds  Stress  transport 
equations  in  the  frame  of  RANS  equations.  The 
treatment  of  incompressibility  (i.e.  the  pressure 
resolution)  appears  in  this  case  more  important  than  the 
respect  of  the  inf-sup  condition. 

2.3.  Analysis  of  numerical  choices  on  LES  of 
channel  flow 

Channel  flow  computations  represent  the  usual  bench¬ 
mark  for  LES  codes.  The  "minimal  flow  unit"  was  first 
introduced  by  Jimenez  etai  (11):  the  computational  box 
is  large  enough  to  sustain  turbulence  but  contains  only 
one  or  two  wall  structures  («  streaks  »).  The  stream  wise, 
wall-normal  and  spanwise  directions  are  noted 
respectively  x,  y  and  z. 


The  computational  domain  is  [2.55h,  2h,  0.9h],  h  being 
the  half  width  of  the  channel.  Periodicity  is  imposed  in 
the  X-  and  z-directions.  The  grid  of  our  LES  is 
composed  by  157  925  nodes,  dispatched  in  25  identical 
planes  in  the  z  direction  (figure  4). 


Fig  4:  unstructured  mesh  in  (x,y)  plane  for  channel  LES. 


The  present  LES  results  are  compared  to  the  DNS  of 
Boudjemadi  (13).  The  Smagorinsky  model  (Cs=0.065) 
with  a  Van  Driest  damping  function  was  used.  Since  the 
ratio  vjv  remains  smaller  than  10%,  this  is  a  quasi 


DNS,  i.  e.  a  good  test-case  for  the  numerical  method. 
Three  numerical  methods  were  tested: 

A:  characteristics  scheme  and  Pl-isoP2  element, 

B:  centred  AB-CN  scheme  on  Pl-isoP2  element, 

C:  centred  on  collocated  PI -PI  element. 

The  spatial  averaging  is  restricted  in  an  unstructured 
mesh  to  narrow  bands  parallel  to  the  wall.  Moreover, 
since  the  computational  box  contains  few  structures:  a 
long-time  integration  is  necessary  to  obtain  reliable 
statistics.  The  time  averaging  was  carried  out  upon  300 
time  units  (LETOTs). 

2.4.  Time  scheme  advancement 

Figure  7  shows  the  time  history  of  maximum  velocity. 
Scheme  B  captures  high  frequencies  oscillations,  which 
are  characteristic  of  turbulent  flow,  whereas  scheme  A 
provides  a  smooth  evolution  due  to  its  low  accuracy  in 
time.  The  rms-value  of  fluctuating  velocities  are  drawn 
on  figure  5.  The  numerical  viscosity  of  scheme  A 
(characteristics)  damps  the  normal  fluctuations,  which 
are  severely  underestimated.  Schemes  B  (centred  +  PI 
isoP2)  and  C  (centred  +  collocated)  seem  equally 
satisfactory,  the  latter  being  however  somewhat  superior 
in  predicting  the  normal  and  span-wise  fluctuations  at 
the  channel  centre. 


Figure  5  rms  values  for  channel  flow. 

The  snapshots  of  vorticity  (not  shown)  and  fluctuation 
iso-contours  (fig.  6)  reveal  for  scheme  A  excessively 
smooth  structures  near  the  wall.  Since  the  normal 
fluctuations  are  underestimated,  the  turbulent  structures 
do  not  migrate  towards  the  centreline  as  well  as  for 
schemes  B  and  C.  Scheme  «  A  »  does  not  capture  the 
essential  features  of  this  turbulent  flow  and  will  not  be 


used  for  subsequent  LES  simulation.  Reasons  for  choice 
of  scheme  C  instead  of  B  is  explained  in  the  following 
section. 
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Fig  6.  Instantaneous  streamwise  and  spanwise 
velocity  contours  in  (y,z)  plane. 


2.5.  Driving  pressure  gradient 


In  LES  of  channel  flow,  it  is  often  felt  more  natural  to 
conserve  the  flowrate  by  driving  the  flow  by  means  of  a 
mean  pressure  gradient,  rather  than  artificially  rescaling 
the  instantaneous  inlet  velocities.  The  mass  flux  is  here 
sustained  by  a  mean  pressure  gradient  At  each  time 
step,  F\  is  function  of  the  variation  of  mass  flux  (13): 


F|"^' 


AtL^h 


with  2  the 


flux,  the  imposed  flux  and  the  length  of  the 

computational  box  in  the  homogenous  direction. 

Fig.7  shows  the  time  history  of  this  external  force.  For 
scheme  A,  the  balance  between  pressure  gradient  and 
friction  exists.  The  centred  scheme  enhances  the 
treatment  of  the  non-linear  term  and  breaks  up  this 
equilibrium.  By  improving  the  pressure  resolution, 
scheme  C  recovers  this  balance.  As  seen  before  on 
Taylor-Green  vortices,  an  identical  resolution  of 
velocity  and  pressure  catches  the  energy  exchanges.  For 
the  minimal  "flow  unit",  the  mean  pressure  gradient  is 
linked  to  the  mass  flux:  the  poor  behaviour  of  scheme  B 
is  attributed  to  the  weak  mass  conservativity  of  the  PI - 
isoP2  element  in  the  LES  context. 


In  the  frame  of  RANSE,  the  velocity  -  pressure  gradient 
correlation  drives  the  energy  exchanges  between 
velocity  component.  This  term  is  here  evaluated  from 
the  instantaneous  fields  and  is  split  in  two  parts: 

The  perturbation  induced  by  the  forcing  term  is 


compared  to:  ( w  )  (  )  . 


For  scheme  B  this  term  is  of  the  same  order  as  the  mean 
resolved  correlation,  thus  pressure  gradient  correction 
used  for  maintaining  the  flowrate  seems  to  introduce  an 
unacceptable  numerical  artefact.  For  scheme  C  this  ratio 
is  less  than  10%. 
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Fig  7:  Driving  pressure  gradient,  Fy  (left), 
and  maximum  velocity  (right)  as  function  of  time. 


Finally,  figure  8  shows  the  balance  of  the  wall  normal 

fluctuation  equation  (usually,  principally  the  w“  budget 
is  shown  since  the  large  magnitude  of  production  masks 
defects  of  the  other  smaller  terms).  Recovering  the 
dissipation  is  not  expected  in  this  LES,  still  graph  e 
indicates  that  with  scheme  A,  small  scales  are  not  well 
resolved,  this  shows  also  on  the  imbalance,  graph  f 
Naturally,  under  estimation  of  v’  also  affects  the 
turbulent  transport,  graph  b. 

The  final  LES  version  (N3S-LES),  has  thus  deviated 
significantly  from  the  original  code.  It  is  based  on  a 
space-&-time-centred,  collocated  FEM  scheme,  with 
linear  basis  functions  on  tetrahedral  elements.  To 
conclude,  the  surprising  finding  here  is  that  methods 
which  are  usually  recommended  in  F.E.  textbooks,  and 


have  indeed  been  successfully  applied  to  RANSE,  are 
not  suitable  for  LES. 


Figure  8  budget  of  the  wall  normal  Reynolds  stress. 

3.  CROSSFLOW  IN  TUBE  BUNDLE 

The  flow  in  a  tube  bundle  is  of  course  of  great  interest  to 
the  power  generation  industry,  not  only  to  study  the 
performance  of  heat  exchangers.  Some  safety  studies 
involve  vibrations  through  fluid-structure  coupling,  or 
large  scale  temperature  fluctuations  eventually  leading 
to  thermal  stripping.  The  geometry  is  simple,  yet  the 
flow  experiences  complex  strains,  making  this  an 
attractive  test  case. 

^-0  x=llinm 


Figure  9  Sub-channel  in  the  tube  bundle  experiment  of 
Simonin  et  al.  Crosses  indicate  measurement  locations. 


3.1.  A  challenge  for  RANSE  simulations 

This  flow  was  considered  for  the  ERCOFTAC/IAHR 
Workshop  on  Turbulence  Modelling  held  at  UMIST, 
June  15-16,  1993.  The  experiment  of  Simonin  & 
Barcouda  (14)  provided  a  database  of  mean  velocities 
and  Reynolds  stresses  for  the  flow  through  an  staggered 
array  of  tubes.  The  tube  diameters  are  D=21.7  nun,  and 
the  distance  between  in-line  cylinders  is  L=45.mm.  The 
Reynolds  number  based  on  the  bulk  velocity  in  the 
subchannel  in  between  tubes,  and  tube  diameter  is 
40000. 

A  variety  of  models  (low  Re  k-e  models,  Reynolds 
Stress  Transport,  both  «  standard  »  and  «  Realisible  » 
were  applied  (appendix).  Some  of  the  conclusions,  as 
presented  by  K.  Hanjalic  &  M.  Leschziner  (15)  are  the 
following  : 

-«  Velocity  profiles  are  predicted  generally  reasonably 
well,  but  stress  profiles  are  poor,  particularly  in  the 
impinging  region  » 

-«None  of  the  results  show  indisputable  superiority, 
though  every  model  refinement  and  upgrading  seems  to 
make  a  positive  contribution  » 

-«  Disappointing  :  Performances  of  RSM,  though  some 
improvements  are  discernible  (stress  behaviour  in  the 
impingement  region  » 

The  following  workshop,  Lisbon  94,  considered  again 
staggered  tube  bundles  with  additional  experimental 
data  form  Me  Grath,  and  Kelemenis,  but  the  former 
provided  only  near  wall  data,  while  some  inconsistancy 
as  regards  mass  conservation  were  found  in  the  latter. 
Again  performances  of  RANSE  models  were 
disappointing,  while  interestingly,  results  using  a  lattice 
Boltzman  scheme  (Eggels)  compared  fairly  well  to  the 
near  wall  data. 

Sebag  (16)  devoted  his  thesis  to  the  application  of 
standard  and  elaborate  second  moment  closures  to  this 
flow,  and  concluded  that  the  major  improvement  in 
predicting  the  Reynolds  stresses  came  from  a  quite 
simple  correction  to  the  dissipation  equation 
(« sensitising  dissipation  to  anisotropy  invariant »  as 
proposed  by  Craft  &  Launder)  and  not  from  the 
modelling  of  the  pressure  strain  term  at  cubic  order. 
Meyer,  providing  his  own  experimental  database  which 
includes  heated  tubes  also  pointed  out  at  the  limitations 
of  second  moment  closures  for  this  tube  bundle  flow, 
attributing  part  of  the  difficulties  to  the  very  high 
turbulence  intensities  in  this  type  of  flow  (35%). 

3.2.  LES  computations 

Finally,  after  the  validation  of  N3S-LES  it  seemed 
natural  to  consider  again  the  experiment  of  Simonin  et 
al). 

The  computational  domain  is  sketched  on  fig.  10.  It 
corresponds  to  the  minimum  domain  for  the 
instantaneous  field  without  enforcing  artificial 
symmetries.  Results  will  show  that  no  coherent  large 
structure  (that  might  have  been  generated  by  periodicity 


on  a  short  domain)  actually  appears.  A  2D  mesh  is 
generated,  structured  in  the  boundary  layers,  and  free 
(Voronoi)  in  between  cylinders.  33  identical  planes  are 
then  extrapolated  in  the  z  direction,  then  the  prisms  are 


The  Reynolds  number  in  the  LES  computation  based  on 
bulk  velocity  in  the  narrowest  section  was  set  to  16  000. 
With  a  mean  friction  velocity  of  u*=  7  mm/s,  the  first 
grid  point  varies  from  y+=1.5  to  4  wall  units  on  the 
circumference  of  the  cylinder.  Because  local 
instantaneous  values  can  be  significantly  higher,  wall 
functions  are  activated  whenever  needed.  The  distortion 
of  the  Reynolds  number  (16  000  instead  of  40  000)  is 
expected  to  have  little  influence  in  this  flow  where  the 
turbulent  intensity  is  35%  (in  opposition  to  flow  around 
a  single  cylinder). 


Quart  du  maillage  plan 


0  X(tnm) 
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Figure  1 1  View  of  1/4^^  of  a  2D  section  of  the  mesh. 


A  first  run  was  computed  with  the  Smagorinsky  model 
(MS  simulation)  and,  after  initialisation,  statistics  were 
accumulated  over  1.745s,  i.e.  three  domain  flow-through 
(this  was  found  sufficient  since  spanwise  averaging  and 
symmetries  are  used  while  gathering  statistics).  For  the 
second  run,  with  the  localised  dynamic  procedure  of 
Piomelli  and  Liu  (16),  statistics  were  computed  over  5 
domain  flow-through  times  (3.468  s). 

The  value  of  the  Smagorinsky  constant  was  set  to  Cs= 
0.065,  the  dynamic  model  later  yielded  somewhat  higher 
values  of  this  coefficient  on  average,  with  of  course 
great  spatial  variability:  high  values  in  the  separated 
shear  layer,  and  slightly  negative  values  very  near  the 
wall  on  the  upstream  half  of  the  tube. 


x(min) 

Figure  12,  Iso-contours  of  the  SGS  viscosity 
obtained  with  the  dynamic  model. 


3.3.  Statistical  results 


For  the  mean  velocity,  the  agreement  of  the  LES  with  all 
available  data  is  excellent,  and  the  results  from  both 
subgrid-scale  models  cannot  be  distinguished.  As 
concerns  the  Reynolds  stresses,  the  agreement  is  also 
good,  so  we  concentrated  on  the  two  section  for  which 
discrepancies  concerning  RANSE  models  were  most 
amplified.  Figure  13  shows  the  wake  to  impingement 
axis  and  the  cross  section  just  behind  the  tube.  The 
striking  and  unusual  feature  on  the  wake  axis  is  that  the 
transverse  velocity  fluctuation  is  far  larger  than  the  axial 
one.  Thus  we  approach  the  stagnation  point  with 

<v^  ,  and  since  -^^<0,  the  productions  of  the 
dx 


stresses  and  kinetic  energy  are  : 


Figure  13  First  and  second  moments, 
solid  lines-LES  with  Smagorinsky  model, 
dashed  lines-LES  with  dynamic  model., 
symbols-experiment  (see  fig  9) 

We  notice  from  the  experiment  and  LES  that  indeed 
is  fairly  constant  (production  balancing  dissipation) 
while  is  decreasing  rapidly  under  the  added  effects 

of  dissipation  and  negative  production.  Since  <v"  , 
even  the  kinetic  energy  is  decreasing,  in  opposition  to 
what  is  found  on  a  stagnation  point  of  a  single  bluff 
body,  or  an  impinging  jet.  Very  near  the  stagnation 

point,  the  cinematic  blockage  of  the  wall  forces  «” 

fluctuation  to  be  converted  into  («  wall  echo  effect  », 
and  the  Dynamic  model  seems  to  better  capture  the  peak 

in  ,  though  shown  only  by  one  experimental  point. 


Figure  14  production  on  the  stagnation  point  axis. 


The  sharp  rise  in  v"  results  in  a  sharp  dip  in  Figure 
14  shows  the  dramatic  overestimation  of  production 
given  by  the  second  moment  closure  (note  the  different 
scaling).  The  experiment  shows  slightly  negative  values, 
but  data  is  not  available  very  close  to  the  wall  to  confirm 
the  sharp  dip  in  shown  by  the  LES. 

In  the  cross  section  just  after  the  tube  (Fig.  13,  bottom) 

at  X=1 1  mm,  the  first  peak  in  v“  is  due  to  the  detaching 
shear  layer.  The  second  peak  corresponds  to  the 
stagnation  point,  and  now  there  are  sufficient 
experimental  points  to  confirm  that  the  trend  given  by 
the  dynamic  model  is  more  accurate  than  that  of  the 
Smagorinsky  model.  On  figure  12  we  had  noted  that  the 
main  difference  between  the  Smagorinsky  and  Dynamic 
model  is  that  the  latter  predicts  zero,  if  not  negative, 
eddy  viscosity  on  the  near  wall  front  portion  of  the 
cylinder.  Vanishing  SGS  dissipation  allows  the  build  up 

of  .  On  the  other  hand,  both  models  overestimate 

.  Since  the  turbulent  viscosity  is  equal  to  five  times 
the  fluid  viscosity  in  most  regions,  some  effect  of  the 
subgrid-scale  model  was  expected,  but  in  any  extent, 
this  is  much  smaller  than  the  extreme  variability 
observed  in  the  aforementioned  workshops  on  RANSE 
models.  From  figure  13,  it  cannot  be  conclude  from  the 
mean  velocity  U,  that  the  recirculation  bubble  is 
accurately  predicted  (for  lack  of  data).  But  further 
downstream,  at  X=  16.5  mm  (figure  15)  we  again  see 
excellent  agreement  with  the  experimental  data  for  both 
models  on  the  U  and  V  components.  Notice  also  the 
very  large  values  of  the  w'  fluctuations  compared  to  u' 
and  v'  ,  especially  on  the  boundary  of  the  upper 
cylinder. 

3.4.  Instantaneous  velocity  fields 

Visualizations  of  instantaneous  iso-streamwise  velocities  (left)  show  that 
large  structures  of  low  velocity  (in  blue/light  grey,  in  circle)  are  detached 
from  the  wake  and  transported  towards  the  stagnation  point  (flow  is  from  left 
to  right).  High  velocity  blobs  (red/dark  grey,  in  ellipse)  occasionally  make 
incursions  in  this  region.  An  explanation  for  the  difficulties  of  standard 
RANSE  models  may  be  that  in  this  wake-to  gap  region,  the  fluctuating 
velocities  are  larger  than  the  mean.  It  is  thus  difficult  to  reconstruct  the  flow 
characteristics,  from  the  sole  mean  velocity.  In  particular,  the  so  called  «  wall 
echo  effect »  introduced  in  second  moment  closures  were  devised  to 
represent  the  «  splatting  »  of  small  eddies  against  an  infinite  (compared  to  the 
integral  scale)  plane  wall.  In  the  present  case,  as  seen  form  the  instantaneous 
fields,  the  size  of  the  structures  is  comparable  to  the  cylinder  size.  For 
instance  the  eddies  are  large  enough  to  «  feel  »  the  curvature  of  the  wall,  i.e. 
at  the  stagnation  point  one  could  imagine  that  the  wall  normal  fluctuations 
are  more  easily  transferred  into  lateral  fluctuations  because  of  this  strong 
curvature. 

Finally,  we  have  pointed  out  the  relatively  large  velocity  fluctuations  in  the 
homogeneity  direction. 

iso-streamwise  velocities. 


X=16.5  mm 


Figure  15  (see  fig.  13) 


Particles  emitted  from  a  vertical  line  (bottom  left)  are 
spread  out  in  the  whole  spanwise  dimension,  showing 
the  fully  three-dimensional  nature  of  the  flow.  This  flow 
is  probably  an  ideal  case  for  the  modelling  the 
mysterious  «  wall  echo  »  terms.  Indeed  they  play  here  a 
dominant  role,  since  overall,  the  w’  fluctuations  (for 
which  there  is  no  direct  production)  are  significantly 
larger  than  the  streamwise  fluctuations  which  are 
alternatively  generated  by  shear  and  damped  by 
cinematic  blocking. 


4.  CONCLUSIONS  AND  PERSPECTIVES 


It  has  been  shown  that  a  standard  « industrial  »  code 
cannot  be  straightforwardly  used  for  LES  which  requires 
specific  dicretizations.  The  low-order  FE  scheme  was 
tested  on  turbulent  flows.  This  analysis  underlines  that 
the  "safest"  element  is  not  adapted  to  LES:  this  statement 
was  suggested  by  the  literature  on  Finite  Element"' 
but  the  limitations  appear  clearly  on  LES  of  turbulent 
flows. 

The  simulation  of  the  flow  in  a  tube  as  a  first  application 
of  LES  in  complex  geometry  yielded  excellent 
agreement  with  experimental  data,  whereas  previous 
RANSE  simulations  had  failed.  This  is  not  «  beginner’s 
luck  »,  LES  was  successful  because  the  size  of  the  larger 
eddies  are  of  similar  order  as  that  of  the  obstacles, 
whereas  RANSE  models  encounter  limitations  precisely 
when  the  integral  length-scale  of  the  turbulence  is 
comparable  with  that  of  the  mean  flow  inhomogeneity 
(this  interpretation  was  previously  proposed  concerning 
the  effect  of  the  Reynolds  number  on  the  recirculation 
behind  a  step  (19)).  Thus  it  may  seem  that  to  much 
effort  has  been  devoted  to  academic  flows,  e.  g  ,  channel 
flows,  where  RANSE  models  yield  excellent  results  for 
tremendously  less  efforts.  Thus  there  seems  to  be  a 
complementary  potential  between  LES  and  RANSE  and 
they  should  not  be  in  considered  in  competition  (just  as 
there  is  no  point  in  insisting  on  using  the  full  3D  Navier 
Stokes  equations  for  coastal  engineering  applications 


instead  of  the  shallow  water  equations  (on  France’s 
shallow  coasts).  Similarly,  LES  of  the  boundary  layer  of 
an  airfoil  at  low  incidence  will  remain  out  of  reach  for 
very  long  (if  not  pointless),  given  the  smallness  of  the 
turbulent  structures  compared  to  the  airfoil,  and  given 
the  reliability  of  RANSE  models  for  such  flows.  On  the 
other  hand  small  bluff  bodies,  e.g.  electronic 
components  in  a  highly  turbulent  channel  flow,  as 
considered  in  the  1997  ERCOFTAC/AIHR  workshop 
does  seem  an  LES  application. 

The  follow  up  of  the  present  simulation  was  logically 
heat  transfer,  but  a  new  project,  in  collaboration  with  the 
French  atomic  energy  commission  (CEA)  has  just 
started  with  the  objective  of  developing  an  LES  code  for 
unstructured  meshes  for  massively  parallel  computers 
(T3E).  The  aim  is  to  perform  LES  of  complex 
geometries  with  over  10  million  nodes,  and  the  tube 
bundle  test  case  will  be  considered  again,  hopefully 
providing  this  time  a  complete  database  (due  begining  of 
1998)  of  the  budgets  of  the  second  moment  closures, 
including  heat  transfer. 
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Appendix  :  RANSE  1993  Workshop  Results  (15) 

(all  graphs  by  T.  Craft) 


A,B,C  :  Various  Low-Re  2  equation  models 
The  3  sets  of  data  (UMIST)  corresponding  best  to  the 
exp  ;  for  the  mean  velocity  (A)  correspondingly  give  the 
highest  turbulent  kinetic  energy  levels. 


D-E :  RNG  models.  Models  were  checked  to  be 
identical,  the  large  differences  observed  are  attributed  to 
a  bifurcation  :  as  the  kinetic  energy  level  is  reduced  by 
the  RNG  correction,  the  flow  is  allowed  to  stagnated  in 
between  tubes,  with  little  mixing  from  the  principle  flow 
passing  the  tubes  with  a  minimum  meandering  effect. 
Less  mean  flow  deformation  in  turn  yields  less 
turbulence  production. 

F,  G,  H  :  Second  moment  Closures 
Standard  IP  model  seems  to  yield  «  best  »  results 

yet  is  severely  underestimated 


:  I 

I,  J  :  cross-section  after  the  tube  (cf.  fig.  15) 
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is  still  severely  underestimated  in  the  wake. 
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ABSTRACT 

The  transient  method  of  measuring  heat  transfer  coefficients  that  uses  liquid  crystals,  since  its 
beginnings  in  the  early  1980s,  has  become  one  of  the  best  ways  of  determining  full  surface 
distributions  of  heat  transfer  coefficient.  The  turbomachinery  research  group  at  Oxford  has 
concentrated  on  the  application  of  the  method  to  numerous  mechanical  engineering  thermal  problems 
specific  to  the  jet  engine.  The  paper  summarises  some  of  the  recent  developments  in  the  technique 
including  the  implementation  of  an  elegant  way  of  producing  a  change  in  the  fluid  temperature.  Recent, 
high-density  heat  transfer  coefficient  measurements  are  discussed  and  a  means  of  integrating  the 
measurements  into  finite  element  software  for  subsequent  data  analysis  is  presented.  The  paper  should 
be  of  interest  to  engineers  interested  in  using  the  most  modern  heat  transfer  measurement  techniques  in 
their  research  and  development  programmes. 


1.  BACKGROUND 

The  transient  method  of  measuring  heat 
transfer  to  thermally  insulating  models  of 
mechanical  engineering  components  has 
become  established  as  one  of  the  most 
effective  ways  of  measuring  complete 
distributions  of  local  heat  transfer  coefficient. 
Its  advantages  include  the  provision  of  data 
over  the  complete  surface  in  one  experiment; 
its  high  accuracy  and  the  ease  of  model 
manufacture.  The  accuracy  depends  on  both 
the  experimental  conditions  and  the 
sophistication  of  the  image  processing  applied 
to  the  liquid  crystal  colour  play.  In  its  most 
straightforward  implementation,  the 
uncertainty  in  heat  transfer  coefficient  is  on  the 
order  of  6%,  see  Baughn  et  al.  (2),  though 
image  processing  that  uses  the  fijll  time 
content  of  the  video  data  can  improve  this. 
When  local  heat  transfer  coefficients  processed 
in  an  image  are  referenced  to  an  accurately 
known  heat  transfer  coefficient,  the  accuracy 
in  heat  transfer  coefficient  ratio  can  be  very 
high  (on  the  order  of  2%)-  see  Wang  et  al.(22). 

The  models  are  made  from  a  good  thermal 
insulator,  normally  perspex,  and  only  require 
the  application  of  a  limited  number  of  surface 
thermocouples  for  crystal  calibration  and 
determination  of  the  initial  temperature.  This 
makes  the  method  particularly  suitable  to  the 
measurement  of  heat  transfer  coefficients  on 
complex  shaped  surfaces.  Chyu  et  al.  (3)  used 
the  method  to  measure  the  heat  transfer 
coefficient  to  dimples  used  to  promote  heat 
transfer  in  internal  cooling  passages.  Wang  et 
al,  (21)  investigated  the  heat  transfer 
coefficient  distribution  over  an  isolated 
pedestal  (or  pin  fin)  with  fillet  radii.  Recent 


examples  of  the  technique  applied  to  external 
heat  transfer  coefficients  include  Robertson  et 
al.(14)  and  (15). 

2.  ESTABLISHMENT  OF  THE 
THERMAL  TRANSIENT 

2.1  Mesh  heater  concept 

Numerous  strategies  have  been  used  to 
establish  the  temperature  difference  between 
the  fluid  and  model  surface.  Early 
experimental  arrangements  included  electrical 
heaters  that  took  time  to  stabilise  at  the 
required  temperature,  (4).  The  heat  transfer 
measuring  method  is  better  conditioned  if  the 
fluid  temperature  undergoes  a  sudden  change 
at  the  start  of  the  test.  For  this  reason,  the 
heater  time  constant  was  overcome  by  using 
plumbing  which  allowed  the  flow  to  by-pass 
the  test  section  whilst  the  former  heated  up. 
Fast  acting  valves  were  then  switched  to  start 
the  heat  transfer  experiment.  Other  approaches 
include  rapidly  inserting  a  pre-heated  or  pre¬ 
cooled  model  into  the  wind  tunnel.  The  large- 
scale  annular  cascade  facility  at  Oxford  uses  a 
pair  of  shutters  to  isolate  the  instrumented 
section  of  the  annulus  from  the  substantial 
flow,  Martinez-Botas  et  al.  (8).  Figure  1,  that 
spans  the  fluid  supply  duct,  Gillespie  et  al  (6). 

2.2  Applications 

The  system  has  been  used  for  research  into 
turbine  blade  cooling  systems,  (5)  and  (20),  for 
investigating  extended  heat  transfer  surfaces 
for  cooling  engine  components  and  for  film 
cooling  experiments  (1).  The  largest  mesh 
heater,  Figure  2,  to  date  consumes  more  than 
40kW  and  is  used  to  heat  a  high  flow  rate 
tunnel  (~10m^/s)  at  speeds  ranging  from  5  to 
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30m/s.  This  work  has  been  reported  in  Neely 
et  al.  (11)  and  led  to  the  optimisation  of 
extended  surface  shapes  used  to  cool  the 
outside  of  engine  components-  see  Robertson 
et  al  (14).  In  this  case,  the  electrical  power 
consumption  was  reduced  by  supplying  current 
to  only  the  central  section  of  the  duct.  The 
powered  section  was  then  surrounded  by  two 
un-powered  meshes  and  analysis  was 
performed  to  ensure  that  the  unheated  flow  did 
not  encroach  on  the  centrally  positioned  heat 
transfer  model.  A  small  (~lmm)  gap 
prevented  current  passing  from  the  central  to 
the  surrounding  meshes.  Figure  3  shows  the 
velocity  measured  340mm  downstream  of  the 
mesh  and  clearly  shows  the  small  jet  caused  by 
the  gap  between  the  unheated  and  the  heated 
meshes.  The  affect  on  the  flow  temperature 
can  be  seen  in  Figure  4  where  the  unheated  and 
heated  zones  are  clearly  denoted.  The  small  jet 
temperature  will  be  between  that  of  the  heated 
and  unheated  zones  and  does  not  show  up 
clearly  in  this  figure. 

Carefril  engineering  of  the  mesh  carrier 
ensures  that  the  temperature  of  the  wires  is 
uniform  so  that  a  uniform  flow  temperature  is 
achieved  downstream.  An  image  representing 
the  mesh  temperature  is  included  as  Figure  6. 
The  mesh  temperature  is  uniform  to  within 
1 .3°C  in  this  case  where  the  flow  temperature 
is  being  increased  by  20°C.  The  mesh 
convective  efficiency,  defined  as 

T  -T 

downstream  upstream 

T  -T 

wire  upstream 

depends  on  the  air  speed,  but  in  this  case  (at 
5m/s)  T]  is  approximately  0.55.  The  non¬ 
uniformity  in  downstream  temperature  can  be 
calculated  as  3.5%  of  the  temperature  increase. 

The  most  recent  meshes  have  been  used  to 
heat  moderately  high  pressure  (-lOOpsi)  air  for 
further  impingement  research,  Figure  7.  For  a 
fixed  mass  flow  rate,  heating  the  air  at  higher 
upstream  density  is  an  effective  means  of 
reducing  the  mesh  face  velocity  (  for  fixed 
area)  and  this  allows  the  heater  size  to  be 
reduced.  The  heater  mesh  in  the  figure  is 
being  used  to  rapidly  switch  the  flow 
temperature  in  a  programme  of  impingement 
cooling  research. 

2.3  Mesh  heater  advantages 
The  mesh  used  has  a  through  flow  to  projected 
area  quotient  of  0.38  and  can,  under  correct 
lighting  conditions,  be  transluscent.  Figure  8 
shows  an  arrangement  used  to  measure  heat 
transfer  coefficients  to  engine  components 
where  the  flow  is  from  a  large  nozzle  through  a 


polycarbonate  perforated  plate(lO).  Figure  9 
shows  an  image  of  the  liquid  crystal  coated 
perspex  model  taken  from  a  miniature  camera 
position  inside  the  nozzle  assembly.  The 
camera  looks  through  the  mesh  and  the  use  of 
a  wide  camera  aperture  effectively  removes  the 
mesh  from  the  image.  The  calculated  heat 
transfer  coefficient  distribution  is  presented  in 
Figure  10.  The  same  approach  has  been  used 
by  Gillespie  et  al.  (6)  to  measure  the  heat 
transfer  coefficients  on  the  upstream  surface  of 
an  impingement  plate,  Figure  1 1. 

All  of  the  electrical  power  supplied  to  the 
mesh  heats  the  air  so  that  an  accurate 
measurement  of  the  average  temperature 
increase  can  readily  be  achieved  by  measuring 
the  power  supplied  to  the  mesh. 


^  downstream  ^  upstream 

mcp 

This  approach  has  proved  very  useful  in  low 
speed  flow  situations  where  the  time  response 
of  the  gas  measuring  thermocouple  is 
prohibitively  slow.  Figure  12  shows  the  level 
of  agreement  between  the  calculated  and 
measured  temperatures  for  recent  blade 
cooling  passage  experiments,  Tsang  and 
Ireland  (18). 

In  some  situations  it  may  be  desirable  to 
generate  a  temperature  profile  using  unheated, 
heated  combinations  of  the  mesh  in  series. 
This  could  allow  profiles  that  have  an 
artificially  thick  cool  region  near  the  wall  to  be 
produced. 

Heat  transfer  measurements  in 
impingement  systems  are  difficult  for  a 
number  of  reasons.  The  local  heat  transfer 
coefficient  changes  by  a  large  factor  from  a 
peak  at  the  stagnation  point  beneath  the  jet  to 
significantly  lower  values  away  from  the  jet  in 
the  channel  between  impingement  and  target 
plate.  One  difficulty  specific  to  the  transient 
method  is  the  problem  of  quickly  establishing 
the  flow  at  a  temperature  different  to  the 
starting  temperature.  In  a  typical  impingement 
array,  the  area  ratio  between  the  plenum 
feeding  the  jets  and  the  total  cross-sectional 
area  of  the  jets  means  that  the  flow  velocity  in 
the  plenum  is  very  low.  This  is  the  region  in 
which  flow  switching  needs  to  be  performed 
and  this  leads  to  problems.  Van  Treuren  et  al. 
(19)  used  a  complicated  by-pass  system  that 
passed  heated  air  into  the  plenum  for 
30minutes  prior  to  the  start  of  the  transient 
experiment.  The  test  section  was  prevented 
from  pre-heating  by  bleeding  room 
temperature  air  through  the  model  into  the 
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plenum.  Three  electrically  operated  fast  acting 
valves  were  synchronised  to  initiate  the 
transient  heat  transfer  experiment.  The  flow 
temperature  switching  in  later  impingement 
heat  transfer  experiments  is  now  achieved 
using  a  planar  heater  mesh  fitted  across  the 
plenum  flow.  The  apparatus  no-longer  requires 
any  fast  acting  valves  since  the  experiments 
start  when  power  is  switched  to  die  mesh. 
Prior  to  this  stage,  room  temperature  air  is 
drawn  through  the  heater  mesh  and  model. 

The  velocity  through  the  heater  mesh  in 
impingement  experiments  is  typically  low  and 
a  step  change  in  the  power  supplied  does  not 
cause  the  flow  temperature  to  change  as  a  step. 
The  normal  surface  temperature  response  of 
the  perspex  substrate  under  a  step  change  in 
fluid  temperature  is  given  by  the  familiar 
equation 


X  erfc{ 


jNL) 

■yjpCk 


When  the  heat  transfer  driving  fluid 
temperature  changes  with  a  time  constant,  x, 
the  equation  for  the  surface  temperature  rise 
becomes,  (6), 
pck 


The  above  equation  is  then  solved  numerically 
for  the  local  heat  transfer  coefficients  from  the 
time,  t,  the  surface  takes  to  reach  the  known 
liquid  crystal  temperature. 

In  situations  when  the  gas  temperature 
change  is  neither  a  step  change  nor  a  first  order 
lag,  the  fluid  temperature  should  be  measured 
and  the  measured  signal  used  in  a  numerical 
procedure  that  sums  the  surface  response  to 
multiple  delayed  steps,  Metzger  and  Larson  (9) 
or  multiple  delayed  ramps,  Saabas  et  al.  (16). 


3.  APPLICATION  TO  IMPINGEMENT 
COOLING  RESEARCH. 

Progress  in  the  hardware  used  at  Oxford  has 
centred  around  a  change  in  the  frame  grabbing 
strategy.  The  card  that  digitises  the  video 
recording  is  fitted  to  a  Pentium  PC  and  the 
ditized  data  stored  directly  in  computer  RAM. 
The  computer  has  195Mbytes  of  ram  fitted  and 
is  capable  of  acquiring  62  full  frames.  In 
practice,  not  all  pixel  information  is  required  to 
resolve  accurately  typical  heat  transfer 
coefficient  results.  In  addition,  not  every 
frame  is  required  for  heat  transfer  processing. 
A  typical  heat  transfer  test  lasts  several 
seconds  and  the  data  is  digitised  over  two  or 
three  sweeps.  Each  frame  is  captured  at  a 
particular  time  from  the  start  of  the  heat 
transfer  experiment.  A  time  signal  is  captured 
in  the  frame  and  used  to  reference  the  intensity 
signal  to  the  elapsed  time. 


4.  IMAGE  PROCESSING- HIGH 
RESOLUTION  HTC  DATA 

Figure  13  shows  the  heat  transfer  coefficient 
distribution  on  the  target  plate  beneath  an  array 
of  impinging  jets.  The  work  is  notable  for  the 
high  data  resolution  that  shows  how  repeatable 
the  heat  transfer  coefficient  signature  between 
rows.  The  exit  from  the  channel  is  from  one 
side  and  spent  flow  from  the  early  rows  forms 
a  cross-flow  on  the  later  jets.  The  research  is 
reported  further  in  (12)  and  (17). 

The  image  processing  strategy  here  was  to 
use  a  combination  of  three  different  narrow 
band  liquid  crystals  that  produce  peaks  in 
intensity  at  three  calibrated  temperatures  at 
every  pixel.  Each  image  includes  more  than 
5000  points  and  the  high  data  resolution  can  be 
seen  to  reveal  interesting  features  about  the 
flow  field.  The  repeatability  of  the  bump  in 
heat  transfer  beneath  each  jet  is  very  noticeable 
even  though  each  jet  is  subjected  to  a  differing 
amount  of  cross-flow.  In  this  arrangement, 
three  different  hole  sizes  have  been  used.  The 
first  three  rows  of  holes  have  the  same 
diameter.  The  fourth  and  fifth  rows  have 
larger  diameter  holes.  The  last  row  has  an 
even  larger  diameter.  One  feature  that  seems 
repeatable  is  the  change  in  heat  transfer 
coefficient  slope  that  exists  between  1.1 
andl  .5d  from  the  hole  axis. 

Earlier  liquid  crystal  experiments  at  Oxford 
by  van  Treuren  tested  an  array  with  a  bigger 
spacing  between  holes.  His  data,  at  a  slightly 
lower  resolution,  exhibited  the  same  shaped 
bumps  beneath  the  jets.  Figure  14  shows 
present  data  normalised  by  the  stagnation  point 
heat  transfer  coefficient.  The  results  from  van 
Treuren  have  also  been  included  and,  despite  a 
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hole  spacing  to  diameter  ratio  more  than  twice 
the  present  value,  there  is  remarkable 
similarity  between  the  heat  transfer  coefficient 
variation.  Van  Treuren's  measurements  were 
made  in  a  transient  heat  transfer  facility  in 
which  the  impingement  plate  was  water  cooled 
to  maintain  a  constant  temperature.  The  liquid 
crystal  data  were  processed  to  obtain  both  the 
heat  transfer  coefficient  and  the  adiabatic  wall 
temperature.  The  latter  was  then  expressed  as  a 
dimensionless  temperature  difference  ratio. 
The  present  data  were  fi*om  a  fully  perspex  rig 
in  which  the  impingement  plate  was  allowed  to 
warm  up  during  the  experiment.  The  heat 
transfer  driving  temperature  is  then  the 
temperature  downstream  of  the  heater  mesh  ( 
the  plenum  temperature)  as  can  be  seen  fi*om 
the  heat  flux-  temperature  plot  included  as 
Figure  15.  The  data  here  could  be  applied  to 
systems  in  wliich  the  plate  temperature  is  the 
same  as  the  plenum  temperature. 

5.  DATA  RECOVERY  FROM  LIQUID 
CRYSTAL  EXPERIMENTS. 

The  subject  of  automatically  recovering  the 
liquid  crystal  colour  change  times  from  video 
recordings  of  experiments  has  been  addressed 
by  Robertson,  (13).  The  software  links  the  co¬ 
ordinates  of  a  numerically  defined  model  into 
pixel  co-ordinates  using  an  analytic  model  of 
the  optical  system.  When  linked  to  the  heat- 
transfer-solving  program,  the  software  can 
automatically  extract  the  measured  heat 
transfer  coefficients  in  real  co-ordinates.  This 
has  proved  invaluable  for  finite  element 
analyses  of  real  engine  components  where  the 
actual  temperature  is  influenced  by  lateral 
conduction  occurring  in  a  metallic  component. 

6.  CONCLUSIONS 

The  paper  summaries  recent  progress  in  what 
has  become  an  established  method  of 
measuring  distributions  of  local  heat  transfer 
coefficient.  Since  its  introduction  three  years 
ago,  the  mesh  heater  has  rapidly  replaced  other 
heaters  as  a  very  convenient  (  and  hence  low 
cost)  way  of  producing  the  thermal  transient. 
The  paper  describes  the  mesh  and  gives  several 
recent  applications  of  its  use.  The  translucent 
nature  of  the  mesh  has  allowed  camera  views 
through  the  heater  and  examples  are  given  of 
the  use  of  this  advantage.  TTie  uniformity  of 
the  mesh  has  been  discussed  and  a  means  of 
reducing  the  total  electrical  power  by  including 
un-powered  sections  presented.  The  mesh  has 
been  shown  to  be  suitable  for  heating 
atmospheric  and  super-atmospheric  flow.  The 
mesh  has  been  used  to  research  the  one  of  the 
most  difficult  cooling  systems-  that  of 
impingement  cooling.  In  this  situation,  the 


change  in  driving  temperature  is  achieved 
easily  and  the  needed  for  complicated 
plumbing  to  pre-heat  the  plenum  avoided. 
Progress  in  the  image  processing  has  increased 
the  data  point  resolution.  A  very  exciting 
computational  method  for  determining  the  co¬ 
ordinates  of  model  surface  at  which  the  crystal 
changes  colour  has  been  introduced. 
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8.  NOMENCLATURE 

Cp  specific  heat  at  constant  pressure  for 

air,  J/kgK 

c  specific  heat  of  perspex,  J/kgK 

h  heat  transfer  coefficient,  W/m^K 

I  current.  Amps 

k  thermal  conductivity  of  perspex, 

W/mK 

m  mass  flow  rate,  kg/s 

t  time,  s 

T  temperature,  °C 

V  voltage,  V 

Greek 

r\  mesh  convective  efficiency 

p  density,  kg/m^ 

X  time  constant,  s 

Subscripts 
g  gas 

O  starting 

s  surface 
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Figure  1  Close  up  of  heater  mesh.  Wire 
diameter  is  40pm. 


MO 


100 

MESH  HEATER 
DETAIL 


•II  dimention*  In  mm 


Figure  2.  Heat  transfer  tunnel  used  to 
investigate  the  thermal  performance  of 
extended  surfaces. 


vertical  dsplacement  from  base  of  test  section  (mm) 


Figure  3  Velocity  profile  measured  340  mm 
downstream  of  the  heater  mesh  showing 
the  edge  of  the  heated  and  unheated 
sections. 


Figure  4.  Vertical  temperature  traverse 
measured  370  mm  downstream  of  the 
heater  mesh  in  extended  surface  research 
tunnel 


Figure  5.  Horizontal  temperature  traverse 
measured  370  mm  downstream  of  the 
heater  mesh  in  extended  surface  research 
tunnel  (130  mm  below  the  cenreline). 


Figure  6  IR  image  of  the  heater  mesh  fitted 
to  the  plenum  of  the  rig  used  to  make  the 
heat  transfer  measurements  to  the 
impinging  jets.  The  extent  of  the  mesh  is 
between  the  white  lines  and  the  surounding 
signals  are  reflections  from  the  rig  intake. 
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Figure  7  Heater  mesh  arrangement  used  to 
heat  high  pressure  air.  Pipe  i.d.  is  50mm. 


Pwmoqircruc 
liquid  oysU 


Figure  8  Arrangement  used  to  measure 
heat  transfer  coefficients  to  engine 
components.  Flow  is  from  large  nozzle 
including  a  perforated  plate. 


Figure  9.  Image  of  engine  mount  ring 
through  the  mesh  and  porous  Perspex  plate 
in  fire  test  rig 


Figure  10.  HTC  distribution  (W/M^/K) 
deduced  from  video  history  of  view  through 
the  mesh  and  porous  Perspex  plate  in  the 
fire-test  rig 
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Figure  11  Nusselt  number  contours  on  the 
front  face  of  a  plate  used  to  produce 
impinging  jets.  The  upper  image  is  from  a 
camera  viewing  through  the  heater  mesh. 


Figure  12  Temperature  increase  calculated 
from  power  supplied  versus  measured  fluid 
temperature  increase  for  a  typical  blade 
cooling  experiment.  Note  the  false  origin- 
temperature  upstream  of  the  mesh  is 
approximately  20°C. 
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Figure  13  DetaUed  heat  transfer  coefficients 
on  the  target  plate  beneath  an  engine 
representative  array  of  impinging  jets.  The 
flow  exits  from  one  side  of  the  channel  in  the 
direction  of  the  arrow. 


Surface  Temperature 


Figure  14  Comparison  of  normalised  heat 
transfer  coefficient  for  impingement  arrays 
with  very  different  spacing. 


Figure  15  Heat  flux  versus  surface 
temperature. 


Figure  16  An  example  of  the  integration  of 
a  finite  element  software  package  with 
video  images  of  a  perspex  mode!  of  an 
engine  component.  The  lower  figure  has 
hidden  lines  removed. 
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'Dirbulent  heat  and  mass  transfer  in  a  uniformly  heated  vertical  tube  with  an  ascending  forced  flow  of  air 

and  a  descending  falling  water  film 


P.  An,  J.  Li  and  J.D.  Jackson 
School  of  Engineering,  University  of  Manchester,  UK 


ESsriiSntal  and  modelling  studies  are  reported  of  combined  heat  and  mass  trai^r  to  air  flowmg  upwards  ® 
eSically  heated  tube  in  the  presence  of  a  thin  Ming  water  film  on  the  inside  surface.  Ejqwnmrats  were  condu^d  to 
examine  L  influence  of  water  injection  temperature  and  water  flow  rate  on  the  tempei^e  distnbution  along  the  heaM 
tube  for  several  values  of  power  input  A  simple  semi-empirical  model  was  developed  to  simulate  the  ocjxnments  md 
enable  the  wall  temperatine  distributions  to  be  computed.  The  results  yielded  by  the  m^el  proved  ^ 
agreement  with  experiment.  TTiey  indicate  that  when  the  water  is  suppUed  at  near  ambient  tempera^,  the  do^t 
Sanism  of  heat  removal  is  convection  by  the  cooling  water.  When  the  water  «  supplied  at 

above,  die  dominant  cooling  mechanism  is  associated  with  evaporation  of  vapour  fi-om  the  fl'ee  surface  of  the  Mmg  water 
film. 


1.  Introduction 

In  the  case  of  passive  cooling,  fluid  is  induced  naturally 
through  a  system  simply  as  a  result  of  heat  transfer 
taifing  place  rather  than  by  the  action  of  a  pump  or  a  fan. 
Due  to  the  low  flow  rates  which  are  typical  of  passively 
cooled  systems  the  convective  heat  transfer  process  can 
be  veiy  different  from  that  of  forced  convection. 
Depending  on  the  thermal  conditions  and  the  geometiy 
of  the  system,  diere  could  be  strong  interactions  l»^een 
free  and  forced  convection.  One  of  the  characteristics  of 
such  mixed  convection  is  that  in  heated  upward  flow  the 
effectiveness  of  heat  transfer  may  be  seriously  impaired 
as  a  result  of  turbulence  production  being  reduced  by 
buoyancy  influences  [1].  In  a  recent  experimental  study 
of  this  problem,  Li  [2]  has  shown  that  such  conditions 
can  certainly  be  encountered  in  a  passively  cooled 
system. 

One  apporach  which  has  been  proposed  to 
overcome  the  problem  of  impaired  heat  transfer,  is  to 
enhance  cooling  by  allowing  a  thin  film  of  water  to  run 
down  the  heated  surfece.  Under  certain  conditions,  a 
considerable  increase  of  heat  removal  can  result  fixjm  the 
evaporation  of  water  at  fte  free  surface  of  the  felling 
water  film.  Heat  and  mass  transfer  into  the  iqjward 
flowing  mixture  of  air  and  vapour  take  place 
simultaneously. 

In  the  present  study,  experiments  have  been 
conducted  to  examine  fee  enhancement  of  cooling  by 
feis  means.  The  system  tested  consisted  of  a  unifon^ 
heated  vertical  tube  wife  a  forced  flow  of  air  ast^ding 
inside  it  at  known  flow  rates  (typical  of  vfeat  might  be 
achieved  wife  naturally  induced  flow).  Modelling  studies 
have  been  undertaken  in  conjunction  wife  fee 
experiments.  A  simple  semi-empirical  model  has  been 
developed  to  simiilate  fee  heat  and  mass  transfer 
processes  involved.  Computational  studies  using  fee 
model  have  enabled  a  better  undersmding  to  be  gained 
of  fee  relative  importance  of  fee  various  mechanisms  of 
heat  transfer  involved. 


2.  The  experimental  investigation 
2.1  Experimental  apparatus 

Figure  1  shows  fee  arrangement  of  fee  test  facility.  The 
test  section  was  a  stainless  steel  tube,  of  internal  diameter 
76  mm,  wall  thickness  1.9  mm  and  length  8  m.  The  tube 
wall  was  resistively  heated  by  passing  electric  current 
through  it.  The  power  was  supplied  by  a  variable  voltage 
AC  supply  system.  The  outside  of  the  tube  was  covered 
wife  a  thick  layer  of  thermal  insulation  and  fee  heat 
generated  within  the  tube  wall  was  removed  by  means 
forced  upward  flow  of  air  within  fee  tube  and  a  Ming 
film  of  water  running  down  fee  inside  surfiice  of  fee 
tube.  The  wall  temperature  distribution  was  measured 
by  thermocouples  attached  to  fee  outside  surfece  of  fee 
tube  at  numerous  axial  locations.  At  fee  top  of  fee 
test  section,  water  could  be  sprayed  uniformly  onto  fee 
inside  surfece  of  fee  tube  from  a  centrally  positioned 
multi-jet  nozzle.  This  could  be  supplied  wife  water  at 
chosen  values  of  flow  rate  and  temperature.  The  water 
was  pre-heated  by  means  of  electrical  immersion  heaters 
in  fee  water  supply  line.  Its  flowrate  was  measured  using 
a  rotameter.  Air  from  a  compressor  was  supplied  to  an 
entry  box  which  enclosed  fee  bellmoufe  intake  at  fee 
bottom  of  fee  tube.  The  flow  rate  of  air  was  monitored 
using  a  carefully  calibrated  flowmeter  installed  in  fee 
supply  line. 

2J,  Experimental  conditions  and  procedure 

The  experimental  conditions  covered  were  as  follows: 

Electrical  power  input:  1,35, 3,00, 4,00  and  5,00  kW 
Air  flow  rate:  0.005, 0.010  and  0.015  kg/s 
Water  injection  rate:  0, 0.013  and  0.027  kg/s 
Water  injection  temperature:  20, 40, 55  and  70  “C 

Air  was  supplied  to  fee  test  section  at  a 
temperature  of  20®C.  The  Reynolds  number  at  inlet 
corresponding  to  fee  three  quoted  air  flow  rates  were 


4600, 9200  and  13800. 

In  a  typical  test  wWi  water  film  cooling,  a 
steady  forced  flow  of  air  was  pun:^)ed  dirough  the  test 
secticm  at  the  required  rate  and  injection  of  water  via  the 
spray  nozzle  was  commenced  Once  a  steady  and 
imiform  water  film  had  been  established,  electrical 
heating  was  sqjplied  to  the  test  section  at  a  constant  rate. 
The  tests  ran  for  a  lengthy  period,  during  which  time  die 
temperature  distribution  along  the  tube  was  monitored  at 
regular  intervals  until  a  steady  state  condition  had  been 
achieved.  Some  tests  were  conducted  widiout  any  water 
being  siqiplied. 

23  Presentation  and  discussion  of  experimental 
results 

Figure  2  shows  tube  wall  temperature  distributions  for 
diree  different  air  flow  rates  (0.015,  0.010  and  0.050 
kg/s)  under  dry  conditions  widi  a  power  input  of  1350 
W.  At  the  lowest  flow  rate  fte  tube  reached  a 
temperature  of  300®C  near  the  top.  This  was  the 
arbitrarily  chosen  operating  temperature  limit  imposed 
to  avoid  damaging  die  test  section.  In  that  test,  the 
conditions  were  such  that  influences  of  buoyancy  were 
strong  and  the  effects  of  non-uniformity  of  fluid 
properties  significant.  Some  effect  of  buoyancy  is 
evident  from  the  sh^  of  the  tube  wall  temperature 
distribution  which  is  irregular.  Figure  3  shows  the 
corresponding  distribution  of  tube  wall  temperature  for 
a  flow  rate  of  0.0 1 5  kg/s  with  a  power  input  of  3000  W. 
This  also  reached  the  operating  temperature  limit  of 
300®C  and  so,  tests  with  the  dry  condition  were  not 
performed  at  lower  air  flow  rates  at  that  power.  The 
results  shown  in  Figure  2  and  3  provide  a  base  against 
which  data  obtained  with  water  film  cooling  can  be 
compared.  Figure  4  shows  such  comparisons  for  an  air 
flow  rate  of  0.015  kg/s  with  a  power  input  of  1350W. 
Results  are  presented  there  for  two  different  water  flow 
rates  (13  g/s  and  27  g/s)  and  four  values  of  inlet  water 
temperature,  (a)  20®C,  (b)  40®C,  (c)  55®C  and  (d)  70®C. 
It  can  be  seen  that  in  each  case  water  film  cooling 
causes  tube  wall  temperature  to  be  greatly  reduced. 
However,  the  behaviour  changes  with  maricedly  increase 
of  inlet  water  temperature. 

From  Figure  4(a),  it  can  be  seen  that  at  the 
lowest  value  of  water  temperature  at  inlet  (20®C),  the 
water  first  increases  in  tenqierature  as  it  descends  and 
dien  its  ten^ierature  begins  to  fell.  As  one  would  expect, 
the  rate  at  which  it  increases  in  temperature  is  greater  in 
the  case  of  the  lower  water  flow  rate  because  the  heat 
input  from  fee  tube  to  fee  water  film  is  essentially  fee 
same  for  both  flow  rates.  The  fell  of  temperature  of  fee 
descending  water  is  brought  about  by  cvjqwration  of 
water  from  fee  surfece  of  fee  film.  The  latent  heat 
needed  is  obtained  at  fee  expense  of  fee  internal  energy 
of  fee  water  left  behind.  The  rate  of  evaporation  is  very 
sensitive  to  fee  temperature  of  fee  water  because 
saturation  vapour  pressure  increases  strongly  wife 
temperature.  Therefore,  evaporation  builds  vp  where  fee 
temperature  of  fee  water  rises  as  it  descends  down  fee 
upper  part  of  fee  tube.  Once  fee  effect  of  heat  removal 
due  to  evaporation  becomes  greater  than  fee  heat  input 


to  fee  water  film  from  fee  tube,  fee  tenq)erature  of  fee 
descending  water  begins  to  fall.  The  changeover  occurs 
earlier  in  fee  case  of  fee  test  wdth  fee  lower  water  flow 
rate  because,  in  feat  case,  fee  water  tenq^erature  is 
higher  due  to  fee  rise  of  ten:q>erature  in  fee  iq>perpartof 
fee  tube  being  greater.  The  water  leaving  fee  bottom  of 
fee  test  section  is  hotter  feat  feat  injected  at  fee  top. 
Most  of  fee  heat  generated  in  fee  tube  wall  is  carried 
away  by  fee  water  leaving  fee  test  section.  The 
air/vapour  mixture,  which  rises  up  fee  tube  increasing  in 
temperature  leaves  at  fee  top  carrying  some  heat  wife  it. 
In  fee  \xpper  part  of  fee  tube  fee  conditions  are  such  that 
some  of  fee  ascending  vapour  is  likely  to  condense  onto 
fee  water  film  releasing  heat  as  it  does.  This  would 
increase  fee  tenq^erature  rise  of  fee  descending  water. 
We  see  that  a  variety  of  heat  transfer  mechanisms  and 
some  rather  complicated  processes  can  be  involved  in 
experiments  wife  water  film  cooling. 

Referring  next  to  Figure  4(b),  water 
temperature  at  inlet  40®C,  we  see  feat  fee  behaviour  is 
rafeer  different  The  temperature  of  fee  descending 
water  changes  much  less.  A  small  net  rise  of  temperature 
occurs  in  fee  case  of  fee  high  water  flow  rate  and  feere 
is  a  small  decrease  in  fee  case  of  fee  lower  water  flow 
rate.  Eviration  of  fee  water  film  is  important  over  fee 
whole  length  of  fee  tube  because  fee  water  temperature 
in  this  case  is  everywhere  ^proximately  40®C.  As  one 
would  expect,  fee  extent  to  which  fee  water  is  cooled  as 
it  fells  is  greater  in  fee  case  of  fee  lower  water  flow  rate 
so  fee  two  distributions  of  tube  wall  temperature  do 
diverge  slightly  in  fee  lower  part  of  fee  test  section 
Generally,  however,  fee  influence  of  water  flow  rate  on 
tube  wall  temperature  distribution  is  quite  small. 

Figure  4(c)  shows  results  for  a  water 
temperature  at  inlet  of  55®C.  The  trend  in  terms  of  fee 
increasing  importance  of  evaporation  wife  increase  of 
inlet  water  temperature  is  continued.  For  both  of  fee 
water  flow  rates,  fee  effect  of  evaporation  is  feat  fee 
descending  water  undeigoes  a  steady  fell  in  temperature. 
Again,  fee  rate  at  which  this  occurs  is  greater  for  fee 
lower  water  flow  rate.  The  air/v^ur  mixture  leaves  fee 
top  of  fee  test  section  not  o^  taking  wife  it  fee 
electrical  energy  generated  within  fee  tube  wall  but  also 
heat  which  has  been  released  by  fee  descending  water 
film. 

Finally,  we  consider  fee  results  shown  in  Figure 
4(d)  for  fee  highest  inlet  water  temperature  (70®C). 
Here,  evaporation  is  of  even  greater  importance.  A  steep 
fell  in  fee  tenq^eratuie  of  fee  descending  water  film  is 
very  evident,  particularly  near  fee  top  of  fee  test  section 
where  fee  temperature  of  fee  water  is  highest.  Again,  fee 
effect  is  greatest  in  fee  case  of  fee  test  wife  fee  lower 
water  flow  rate. 

Figure  5  shows  a  corresponding  set  of  results 
for  a  lower  air  flow  rate  (0.01  kg/s)  and  fee  same  power 
input  (1350  W).  The  general  picture  is  very  simOar, 
particularly  for  fee  lower  values  of  water  ten:^)erature  at 
inlet  (20®C  and  40®C).  For  water  temperature  of  55®C 
and  70®C,  fee  tube  wall  tenq)eratures  are  slightly  higher 
than  fee  corresponding  ones  shown  in  Figure  4,  a 
consequence  of  fee  effectiveness  of  heat  and  mass 
transfer  into  fee  air/v2qx)ur  being  less  wife  fee  lower  air 


flow  rate.  However,  Ae  effect  is  quite  small  and  main 
result  of  reducing  Ae  air  flow  rate  is  Aat  Ae  influence 
of  water  flow  rate  is  reduced  for  values  of  inlet  water 
temperature  of  40®C  and  above.  The  results  shown  m 
Figure  5  generally  follow  a  similar  pattern  to  Aose  m 
Figure  4. 

Figures  6,7,8  and  9  show  Ae  influence  of  air 
flow  rate  on  tube  wall  temperature  distribution  for  inlet 
water  temperatures  of  20®C,  40®C,  55®C  and  70®C.  whh 
a  power  mput  of  1350W.  Three  values  of  air  flow  rate 
are  covered  in  each  case. 

For  a  water  temperature  at  inlet  of  20®C  (Figure 
6),  it  can  be  seen  Aat  air  flow  rate  only  affects  Ae  heat 
transfer  process  m  Ae  lower  part  of  Ae  tube.  There, 
evaporation  from  Ae  descendmg  water  film  is 
important.  The  tube  wall  temperature  needs  to  be  higher 
in  Ae  case  of  Ae  lower  air  flow  rate  to  transfer  heat  by 
turbulent  convection  from  Ae  water  film  into  he 
air/vapour  mixture.  The  effect  is  most  apparent  at  Ae 
lowest  air  flow  rate  of  0.005  kg/s  where  some  buoyancy 
induced  impairment  of  turbulent  diffusion  is  beginning 
to  be  experienced.  This  is  of  importance  in  Ae  case  of 
Ae  corresponding  results  for  inlet  water  temperatures  of 
40®C,  55®C  and  70®C,  shown  in  Figure  7,8  and  9.  The 
heat  transfer  mechanisms  present  in  Ae  cases  shown  in 
Figures  6  to  9  are  essentially  Ae  same  as  Aose  discussed 
earlier  wiA  Ae  added  influence  of  buoyancy  at  Ae 
lowest  air  flow  rate. 

We  next  consider  Ae  effect  of  increasing  Ae 
power  input  to  Ae  test  section.  Figures,  9, 10  11  and  12 
togeAer  provide  a  picture  of  what  happens  as  power 
input  is  increased  in  Ae  range  of  1350  W  to  5000  W  for 
a  fixed  value  of  inlet  water  temperature  of  70®C.  Three 
air  flow  rates  are  covered  in  each  case  (0.005  kg/s, 
0.010  kg/s  and  0.015  kg/s).  The  most  striking  feature  of 
this  series  of  results  is  that  Ae  tube  wall  temperature 
only  increases  by  a  raAer  modest  amount  as  a 
consequence  of  Ae  power  being  increased  by  a  Actor  of 
over  Aree.  The  peak  wall  temperature  wiA  a  power 
input  of  5000  W  is  only  82®C,  as  compared  wiA  72®C 
for  a  power  input  of  1350  W.  In  each  ease  Ae  dominant 
mech^sm  for  heat  removal  from  Ae  tube  is 
evaporation  of  water  from  Ae  descending  water  film. 
The  turbulent  convection  process  is  raAer  dependant  on 
air  flow  rate  and  a  systematic  increase  of  wall 
temperature  is  evident  in  each  case  as  this  is  reduced.  In 
ail  cases,  not  only  is  heat  removed  from  Ae  tube  but  also 
Ae  falling  water  is  cooled  down  by  a  significant  amount. 
The  system  is,  in  effect,  operating  as  an  evaporative 
cooler.  In  Figure  13  Ae  results  for  a  power  input  3000 
W  are  presented  again.  However,  Ais  time  Aey  are 
plotted  on  a  different  scale  along  wiA  Ae  tube  wall 
ten^)erature  distribution  for  dry  conditions.  This  figure 
highlights  Ae  effectiveness  of  evaporation  from  a  Ming 
water  film  as  a  mechanism  for  heat.  It  is  apparent  Aat 
wiAout  water  film  cooling  Ae  tube  waD  temperature 
would  have  been  greatly  in  excess  of  300®C  for  power 
inputs  of  4000  W  and  5000  W  and  serious  overheating 
of  Ae  test  section  would  have  occurred.  However,  using 
water  film  cooling  wiA  water  at  70®C,  Ae  tube  is 
maintained  at  a  perfectly  acceptable  temperature  which 
is  actually  below  Aat  of  Ae  cooling  water. 


3.  The  modelling  study 

3.1  Physical  system 

In  Ae  modelling  study,  it  was  assumed  Aat  Ae  tube  was 
uniformly  heated  and  perfectly  insulated  on  Ae  outside. 
A  circumferentially  uniform  film  of  water  flowed 
downwards  on  Ae  inside  surface  of  Ae  tube  under  Ae 
action  of  gravity.  Water  at  temperature  Ti,o  is  supplied  at 
Ae  top  of  Ae  tube  at  a  flow  rate  of  mf.  Air  is  supplied  to 
Ae  tube  at  Ae  bottom  wiA  a  mass  flow  rate  m&  at 
ten^)erature  Tg,o  wiA  relative  humiAty  Rh. 

3.2  The  modelling  approach 

The  test  section  is  Avided  into  N  elements  along  its 
lengA  (see  Figure  14).  Each  element  is  made  up  of  three 
regions,  the  wall,  liquid  film  and  gas  flow  regions.  The 
variables  for  each  element  (see  Figure  15)  are  Ae  tube 
wall  temperature  Tw4,  Ae  gas  temperature  Tg^,  Ae 
concentration  of  vapour  in  Ae  gas  stream  dv,i  and  Ae  gas 
velocity  Vg,i,  The  input  parameters  are  Ae  power  input 
Q,  Ae  initial  water  film  temperature  (i.e.  water  injection 
temperature)  Ti,o  (=Tw,n),  Ae  water  injection  flow  rate 
mfi,  Ae  air  flow  rate  m&,  Ae  air  inlet  temperature  Tg.o 
and  air  inlet  relative  humidity  Rh.  The  mass  and  energy 
balance  equations  involving  Aese  variables  and 
parameters  are  formulated  for  each  of  Ae  three  regions. 
The  rate  of  convective  heat  transfer  from  Ae  water  film 
to  Ae  gas  stream  at  Ae  water-air  interface  is  calculated 
using  a  well-established  empirical  correlation  [3]  wiA 
account  being  taken  of  Ae  effect  of  buoyancy  iidluences 
[2].  The  rate  of  mass  transfer  is  calculated  by  making  use 
of  Ae  heat/mass  transfer  analogy.  Then,  wiA  Ae 
boundary  conditions  specified,  Ae  tube  wall 
temperature,  Ae  bulk  temperature  of  gas  and  oAer 
parameters  can  be  computed. 

33  Modelling  results 

Figures  16  to  19  show  Ae  temperature  distributions  for  a 
fixed  water  injection  rate  of  0.027  kg/s  and  a  fixed  water 
temperature  of  70®C  for  four  power  inputs  (1350,  3000, 
4000  and  5000  W). 

In  each  figure,  Ae  Aree  sets  of  experimental 
data  (denoted  by  markers)  correspond  to  mass  flow  rates 
of  air  through  Ae  system  of  0.005, 0.010  and  0.015  kg/s. 
The  corresponding  wall  temperature  distributions 
obtained  using  Ae  Ae  semi-empirical  model  are  shown 
by  M  and  broken  lines. 

3.4  Discussion  of  the  modelling  results 

It  can  be  seen  Aat  Ae  semi-empirical  model  reproduces 
Ae  e?q)erimentally  observed  behaviour  extremely  well 
for  Ae  cases  shown.  Quite  good  agreement  was  also 
obtained  in  Ae  case  of  experiments  wiA  water 
temperature  at  inlet  of  55 ®C.  Agreement  between  Ae 
preActed  results  and  experiment  for  water  temperatures 
at  inlet  of  40®C  and  2b®C  is  less  sadstactory  but 
neverAeless  Ae  main  trends  exhibited  by  Ae 


experimental  data  are  still  reproduced  quite  well. 

There  are  diree  mechanisms  involved  in  the 
cooling  process,  namely,  laminar  convection  in  the 
descending  water  film,  turbulent  convection  of  heat  in 
the  upward flow  of  air  and  vapour  and  evaporation  of  the 
water  film.  The  numerical  simulations  enable  us  to 
evaluate  the  separate  contributions  to  overall  heat 
removal.  Figure  20  shows  a  breakdown  of  heat  transfer 
to  the  water  to  the  air/v2qx)ur  mixture  for  e;q)eriments 
with  a  power  inpvX  of  1350  W  for  the  various  values  of 
water  temperature  at  inlet.  Figure  20(a)  shows  the  rate  of 
heat  removal  due  to  laminar  convection  by  die  water. 
With  low  water  injection  temperatures  (20  and  40®C), 
die  removal  rate  is  positive,  i.e.  the  felling  water  film 
carries  heat  away  widi  it  With  higher  water  injection 
temperatures  (55  and  70®C),  the  value  is  negative.  The 
descending  water  film  leaves  the  test  section  colder  than 
when  it  was  injected.  Figure  20(b)  shows  separately  die 
rates  of  heat  removal  due  to  turbulent  heat  convection  in 
the  mixture  and  evjqxiration.  Generally,  both  get  bigger 
with  increase  in  water  injection  temperature.  It  can  be 
seen  that,  die  rate  of  heat  removal  due  to  evaporation  is 
much  higher  than  that  due  to  turbulent  convection  in  the 
mixture.  This  confirms  that  evaporation  plays  a  dominant 
role  in  die  overall  cooling  process  for  such  conditions. 
From  Figures  20(a)  and  20^),  we  can  see  in  detail  how 
the  contributions  of  these  heat  transfer  mechanisms  vaiy 
with  the  water  injection  temperature.  For  example, 
consider  an  air  flow  rate  of  0.005  kg/s,  a  water  flow  rate 
of  27  g/s  and  a  power  input  of  1350  W.  The  percentage 
contributions  of  the  individual  mechanisms  to  the  overall 
heat  removal  are  as  tabulated  below. 


Table  1  Percentage  heat  removal  by  laminar  convection 
in  the  falling  water  film,  turbulent  convection  of  heat  in 
the  air/vapour  mixture  and  evaporation 


Water 

injection 

temperature 

Laminar 
convection 
by  falling 
water  film 

Turbulent 
convection 
of  heat  in 
mixture 

Ev^oration 
from  the 
water  film 

IQ^C 

90% 

1% 

9% 

40»C 

43% 

9% 

48% 

55®C 

-5% 

12% 

93% 

70»C 

-117% 

18% 

199% 

We  see  diat  with  die  lowest  water  injection 
temperature  of  20®C,  heat  removal  by  the  flowing  water 
represents  90  percent  of  the  overall  heat  removal. 
Evaporation  contributes  only  9  percent  This  pattern 
changes,  however,  as  the  water  injection  temperature  is 
increased.  For  die  highest  water  injection  temperature  of 
70®C,  the  felling  water  film  does  not  remove  any  heat,  in 
feet  it  releases  a  considerable  amount.  A  lot  of  heat  is 
removed  as  a  result  of  evaporation.  The  percentage  of 
heat  removal  by  tinbulent  convection  of  heat  is  much 

lower  and  varies  from  1%  to  18%. 

Figures  21(a)  and  21(b)  show  the  rate  of  heat 


removal  by  the  three  mechanisms  as  a  function  of  power 
input  and  air  flow  rate  for  a  water  injection  rate  of  0.027 
kg/s  and  a  water  tempCTature  of  70®C  at  inlet  It  can  be 
seen  diat  evaporation  is  in  all  cases  die  dominant 
mechanism  of  heat  removal. 

The  modelling  study  reported  hear  has  enabled 
a  very  useful  picture  to  be  obtained  of  the  extent  to 
vriiich  the  various  mechanisms  of  heat  transfer 
contribute  to  the  removal  of  heat. 

4.  Conclusions 

The  ^plication  of  water  film  cooling  substantially 
reduces  tube  wall  temperature  as  compared  to  that  for  air 
flow  alone.  Three  mechanisms  are  involved,  namely, 
l^yninar  convection  in  the  descending  water  film, 
turbulent  convection  of  heat  in  the  upward  flow  of  air 
and  evsqxiration  from  the  water  film  into  die  air/vapour 
mixture.  The  e?q)erimental  and  modelling  studies 
reported  here  have  shown  that  with  water  injection  at  the 
higher  temperatures  used  in  die  tests  (55®C  and  70®C) 
evaporative  heat  transfer  is  the  dominant  heat  removal 
mechanism.  With  low  power  levels  and  low  inlet  water 
temperature,  the  principal  cooling  mechanism  is  laminar 
convection  in  the  descending  water  film. 
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Power  Input:  1 350  W 
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Dietanoe  from  the  bottom  of  the  tube  (m) 


Figure  2  Tube  well  temperature  dietrlbutlon  under  dry  condttione 


Figure  3  TUbe  wall  temperature  dietrlbutione  under  dry  eondittone 
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Power  input:  1350W 

Air  flow  rate:  0.01 5  kg/e,  Inlet  temp:  20  deg  C 


Power  Input:  1350W 

Air  flow  rate:  0.01 5  kg/e,  Inlet  temp:  20  deg  C 


Flgur#  4  TidM  wall  tamparatura  dlatributiona  with  watar  film  cooling 
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Figurt  6  TUb«  wall  tamparatura  dlatrlbutlona  with  watar  film  cooling  Flgura  7  TUba  wall  tamparatura  dlatrlbutlona  with  watar  film  cooling 


Powar  Input:  1360  W 

Water  flow  rmte:27  g/e,  Water  temp;  56  deg  C 
Ambient  temp;  20  deg  C 


Power  Input:  1350  W 

Water  ftow  rate;  27  g/e,  Water  temp:  70  deg  C 
Air  inlet  temp:  20  deg  C 


Figure  8  Ttiba  wall  tampartf  ura  dlatrlbutlona  with  watar  film  cooling  Flgura  8  'Hiba  wait  tamparatura  dlatrlbutlona  with  watar  Aim  cooling 
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Pigurs  12  TUb#  wall  tsmpsraturs  distributions  with  watsr  film  cooling 


PIgurs  13  Tubs  wall  tsmpsraturs  distributions 
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Flgur«14  Moddling  details 
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FIgura  15  Paramatart  and  varfablaa  for  alamant  1-1, 1  and  1+1 


Powar  input;  1350  W 

Watar  flow  rMa:  27  g/i,  Watar  tamp:  70  dag  C 
Air  iniat  tamp:  20  dag  C 


Powar  Input:  3000  W 

Watar  flow  rata:  27  g^,  Water  temp:  70  deg  C 
Ambient  tamp:  20  deg  C 


FIgura  16  Comparison  of  prsdictsd  tub#  wall  t#mp#ratur# 
distributions  with  sxpsrimsnUt  rssults 


Figure  17  Cempsriaon  of  pradletad  tuba  wall  tampsrsturs 
distributions  wtth  sxpartmsntsl  results 
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(b)  Rataa  of  haat  ramoval  by  convactlon  and  avaporatlon 
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